%% process_plot_scalar
function process_plot_scalar(mesh, fespace, SOL, varargin)
x = fespace.node(:,1);
y = fespace.node(:,2);
z = fespace.node(:,3);
if isempty(varargin)
    u = full(SOL);
    figure("WindowStyle", "docked");
    tri = convert_face(mesh, fespace);
    patch("Vertices", [x, y, z], ...
        "Faces", tri, ...
        "FaceVertexCData", u, ...
        "FaceColor", "interp", ...
        "EdgeColor", "black", ...
        "FaceAlpha", 0.7);
    colorbar;
    xlabel("x");
    ylabel("y");
    zlabel("z");
    title("numerical");
    axis("equal", "tight");
    view(120, 30);
else
    figure("WindowStyle", "docked");
    scatter3(x, y, z, 20, SOL, "filled");
    colorbar;
    xlabel("x");
    ylabel("y");
    zlabel("z");
    title("numerical");
    axis("equal", "tight");
    view(120, 30);
end
end

%% convert_face
function new_tri = convert_face(mesh, fespace)
if strcmp(fespace.basis_type, "P1")
    new_tri = mesh.tri(:,2:4);
elseif strcmp(fespace.basis_type, "P2")
    t = [1 5 6 7;
        1 5 6 8;
        1 5 7 9;
        1 6 7 10;
        2 6 10 7;
        3 8 6 5;
        4 5 9 7];
    elem = reshape(fespace.elem(:,t), [], 4);
    Ne = size(elem,1);
    tri_all = zeros(Ne*4,3);
    t = [4 3 2;
        1 3 4;
        4 2 1;
        1 2 3];
    k = 1;
    for e = 1:Ne
        tri_all(k:k+3,:) = reshape(elem(e,t), 4, 3);
        k = k+4;
    end
    [~, ib, ~] = unique(sort(tri_all,2), "rows", "stable");
    new_tri = tri_all(ib,:);
end
end